Quantum search with interacting Bose-Einstein condensates 
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One approach to the development of quantum search algorithms is the quantum walk. A spatial 
search can be effected by the continuous-time evolution of a single quantum particle on a lattice 
or graph containing a marked site. In many physical implementations, however, one might expect 
to have multiple particles. In interacting bosonic systems at zero temperature, the dynamics is 
well-described by a discrete nonlinear Schrodfnger equation. We investigate the role of nonlinearity 
in determining the efficiency of the spatial search algorithm within the quantum walk model, for the 
complete graph. The analytical calculations reveal that the nonlinear search time scales with size of 
the search space N like \/~N, equivalent to the linear case though with a different overall constant. 
The results indicate that Bose-Einstein condensates in optical lattices could be natural systems for 
the implementation of the quantum search algorithm. 
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I. INTRODUCTION 

Searching is arguably one of the most important prob- 
lems in computer science. The search problem consists of 
finding a particular (marked) element in a set containing 
N items. The simplest classical approach is to uniformly 
sample the set until the marked element is found, which 
on average occurs in a time t = O(N). Though the sam- 
pling can be recast into the framework of Markov chains, 
the O(N) scaling of the classical search time is optimal. 

Grover [l|, Q constructed an algorithm based on quan- 
tum mechanics that can find the marked element of a 
set quadratically faster, in a time t = 0(y/~N). This 
was proven to be optimal |3j. The quantum search al- 
gorithm can be recast in terms of quantum walks, the 
quantum mechanical extension to Markov chains. The 
origins of the discrete-time quantum walk can be traced 
to Meyer 0, Q, and the concept was further developed 
in 2001 by several others [6-8]. The continuous-time 
quantum walk, the quantum analog of a continuous- 
time Markov process, was introduced by Farhi and Gut- 
mann [9] and extended by them and Childs [10]. Nu- 
merous algorithms based on the quantum walks were 
soon developed that were shown to be more efficient 
than the best- known classical algorithms [lT| - [l5j . Among 
these are quantum walk search algorithms, based both 
on the discrete-time quantum walk Il6l4l9j and on the 
continuous-time quantum walk [2CH22j . In these ap- 
proaches, the set corresponds to the vertices of a graph, 
and the marked element is one distinguishable vertex. 
While the spatial search time attains the optimal scaling 
on most graphs, it remains at best t = 0(y/N log N) for 
quantum walks on the two-dimensional square lattice de- 
spite much effort QJ, [H H, [H[, and t = 0(N) in one 
dimension. 

Ultracold atoms are promising physical systems for the 
implementation of the quantum walk search algorithm. 
For example, three steps of a discrete-time quantum walk 
were realized with 25 Mg+ ions in a linear Paul trap [25] ; 
longer quantum walks were effected more recently with 
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Quantum walks have also been imple- 



mented with neutral atoms confined in optical lattices 
using Bose-Einstein condensates of 87 Rb 27] as well as 
single 133 Cs atoms [28j. Because the behavior of quan- 
tum walks is governed by quantum interference, it is not 
necessary to restrict physical systems to single atoms: 
Bose-Einstein condensates (BECs) consisting of millions 
of atoms can maintain their phase coherence over long 
timescales [29j, and are simpler to prepare and manip- 
ulate than single atoms or ions. The phase coherence 
times in BECs are limited by particle interactions, but 
these effects can be minimized by reducing the interac- 
tion strength through the use of Feshbach resonances 1301 . 
Moreover, an individual site of an optical lattice [3l| 
could in principle be 'marked' using the recently demon- 
strated single-site addressing [32l ]. 

The implementation of a continuous-time quantum 
walk using a BEC of ultracold atoms in an optical lat- 
tice is equivalent to allowing the bosons to evolve under 
their governing Hamiltonian. Because of particle interac- 
tions, the resulting 'Gross-Pitaevskii' equation of motion 
in the mean-field approximation (which accurately cap- 
tures the dynamics at low temperatures and large parti- 
cle numbers) is nonlinear [33| . In principle, the presence 
of nonlinearity in quantum dynamics can radically alter 
the performance of quantum algorithms 3J] , even allow- 
ing NP-complete problems to be solved in polynomial 
time. One might therefore conjecture that the timescale 
for the quantum search problem could be modified by the 
presence of nonlinearity. That said, the nonlinearity that 
appears in the Gross-Pitaevskii equation has its origins 
in ordinary linear quantum mechanics, which would ap- 
pear to rule out any improvement in the quantum search 
time (though it could always be worse). In any case, it 
is important to know if the nonlinearity in the Gross- 
Pitaevskii equation would preclude the use of BECs for 
the implementation of a quantum search. The results 
presented here indicate that interacting BECs can indeed 
implement a quantum spatial search algorithm. 

The influence of (a physically motivated) nonlinearity 
on the time to effect a quantum walk spatial search is the 
central question addressed in this work. We consider the 
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quantum search algorithm on the complete graph using 
a continuous-time quantum walk based on the discrete 
nonlinear Schrodinger equation. Though the complete 
graph is perhaps difficult to realize experimentally, its 
study offers several theoretical advantages. The search 
time obtained in the linear quantum walk algorithm has 
previously been shown to be optimal [20j . The neighbor- 
hood of every vertex corresponds to all other vertices, so 
that the quantum walk naturally mimics a uniform sam- 
pling of the set elements. Last, the symmetry of the com- 
plete graph allows the vertex set to be decomposed into 
two inequivalent elements: the marked vertex, and the 
set of unmarked vertices. This allows the iV-dimensional 
Hilbert space to be reduced to two dimensions, greatly 
simplifying the analysis of the nonlinear problem. 

This manuscript is organized as follows. In Sec. Hi] 
we review the continuous-time quantum walk approach 
to the spatial search problem, and derive the associated 
nonlinear equation of motion for a quantum search based 
on interacting Bose-Einstein condensates. The analytical 
results are presented in Sec. IIII1 and the criteria for a 
complete search (unit output probability on the marked 
vertex in the limit of large N) and an incomplete search 
are derived. The results are summarized in Sec. IIVI 

II. BACKGROUND 

A. Continuous-time quantum walk search 
algorithm 

In the continuous-time quantum walk, the state of the 
quantum walker \ip) is evolved in time by the action of 
the Hamiltonian Ho = —jL, where L is the Laplacian 
operator and 7 is the transition amplitude. Given an N- 
dimensional graph G — (V, E) where V = {1,2,..., N} 
and E correspond to the set of vertices and edges re- 
spectively, one can define the Laplacian as L = A — D, 
where A = A(G) is the adjacency matrix associated to 
the graph G and D is a diagonal matrix whose elements 
are Da — ^T. Aij — deg(i), the degree of vertex i (the 
inclusion of the diagonal is not needed if the graph is 
regular). The adjacency matrix specifies the graph con- 
nectivity and its matrix elements are defined as 

A^i 1 ^ )eE (1) 
I otherwise. 

In the continuous-time quantum walk, one associates 
each vertex i to a basis vector the set of basis vectors 
spans the 7V-dimesional Hilbert space. The state of the 
quantum walker is \ip(t)) = CLi{t)\i) , where cii(t) are 
time-dependent complex coefficients. The quantum walk 
is then effected by performing the unitary transformation 
\ip(t)} = e trf ot \ip(Q)) on the particle state vector (H is set 
to unity in this work). 

In the continuous-time quantum walk search algorithm 
of Childs and Goldstone [20| . one of the basis vectors \w) 



is treated differently. This is accomplished by introduc- 
ing a marking or oracle Hamiltonian: 

H w = -\w)(w\. (2) 

The quantum state is then evolved according to the total 
Hamiltonian H = H + H w , i.e. \tp(t)) = e l ~< m \ip(0)). It 
was shown that if the initial state is chosen to be the 
uniform superposition of all sites 

1 N 

hH0)) = |S)^-=5>), (3) 

v i— 1 

then there exists a time t s and value of 7 for which the 
probability on the marked site \(w\ip(t s ))\ 2 = \ip w (ts)\ 2 
will attain unity. For the complete graph, t s = ^ v^/V and 
7 = -h . For the hypercube and an m-dimensional square 
lattice for m > 4, the search time remains t s = 0(y/N). 

B. Discrete nonlinear Schrodinger equation 

A convenient starting point for the description of BECs 
in lattices is the Bose-Hubbard Hamiltonian [35| : 

Hbh = -7 E ( h Y b * + + I E h ' fa - !) > ( 4 ) 

where bi (SJ) annihilates (creates) a boson in site i, and 
hi = b\bi is the number operator; 7 > denotes the 
amplitude for a boson to hop between sites i and j, and 
U > is a repulsive on-site interaction energy. The ba- 
sis states for M bosons can be conveniently represented 
in the Fock (number) representation \ni, ri2, . . . , njv) = 

Eli (M^j 1 0) , where |0) is the particle vacuum state and 

n i = M. The Heisenberg equation of motion for the 
bosonic field operators is 

ib k = [b k ,H BE ] = -7 E h + Ufi k b k - —h, (5) 

(i,k)eE 

where [A, B] denotes the commutator of A and B. The 
last (constant) term can not affect the dynamics, so it 
will be ignored. 

At zero temperature when 7 3> U and the number 
of particles M is large, the bosons will form a BEC. 
Mean-field theory then provides an excellent description 
of the particle dynamics. This is equivalent to ignoring 
the quantum fluctuations of the bosonic field operators: 
Sbi = bi — (bi) w 0. Defining the BEC wavefunction 
"fi = (bi), the equation of motion becomes 

% ^fk = -l E ^ t +U\^ k \ 2 ^k- (6) 

(i.k)eE 

Because all of the bosons are now assumed to be in the 
same single-particle state, we have the normalization con- 
dition J2i \^i\ 2 = M- It is convenient to define the BEC 
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wavefunction in terms of a new wavefunction 'I'.; = v Mipi 
so that IV'il 2 = !• Note that the single-particle basis 
states can now be chosen to be \i), so that the BEC wave- 
function is rpi — (i\tp). Defining g = UM, one obtains the 
discrete nonlinear Schrodinger or Gross-Pitaevskii (GP) 
equation: 



d 



(7) 



For graphs with constant degree d the Laplacian is L — 
A—D = A—dl, so the dynamics of the BEC wavefunction 
is equivalently described by the Gross-Pitaevskii Hamil- 
tonian 



N 



ff G p--7i + 5^W 2 |fc)(fc|. 



(8) 



fe=i 



The non-linear quantum walk search Hamiltonian then 
takes the following form: 



H = H, 



GP- 



\w)(w\ 



N 



- 7 L-\w){w\+g^\^\i){i\. (9) 



III. NONLINEAR QUANTUM WALK SEARCH 
ON THE COMPLETE GRAPH 

A. Reduction to two dimensions 

Consider the action of the Hamiltonian (|9]) for the N- 
vertex complete graph . The complete graph is associ- 
ated to the adjacency matrix with elements Atj = 1 — 8%j. 
The vertex degree is constant d = N — 1 so the Laplacian 
is L = A - (JV - 1)1 = J - NI, where J is the JV- 
dimensional all-one matrix. In terms of the initial state 
\S) defined inEq. ©, the Laplacian is L = N\S)(S\-NI. 
Since —NI corresponds to a constant energy shift it can 
not change the dynamics of the system. The Hamiltonian 
for the quantum walk search @ then becomes 



JV 



h = - 7 n\s){s\ - \w)( w \ m 2 m\ 



(10) 



1=1 



In the absence of the nonlinear term, this Hamiltonian 
corresponds to a two-level operator for the states \S) and 
\w). The Hamiltonian then rotates \S) into \w) in time 
t s ~ VN inversely proportional to the two states' overlap 
(S\w)=l/VN. 

Though the Hamiltonian ()10j) is nonlinear, it would be 
desirable if it could also written as a two-level operator. 
For the complete graph, the initial condition ^ is 



Mo)) = |s> = -^y>> = J=i 

Vn^ Vn 



N- 1 



N 



a), (11) 



where \a) labels any state orthogonal to the marked state 
\w). The initial amplitude is the same on any site Ice), and 



the nonlinear part of the search Hamiltonian preserves 
this symmetry: 

£><| 3 |i)<i| = \Tp w \ 2 H(w\ + (N-l)\i> a \ 2 \a)(a\. (12) 



The search Hamiltonian, including the nonlinear part, 
therefore remains a two-level operator on the states 
spanned by \w) and \a). 

The equations of motion can be found using 

{w\i-\iP) = {wm); (a\i-\i>) = {a\H\i>). (13) 

These yield 

d 

l 7T^™ = -i{N - l)ip a - (1 + j)ip w + g\ip w \ ipw; (14a) 
at 



1 1T^<* = ~^( N ~ - l^w + g\ip a \ 2 ip a - (14b) 
at 



The dynamical equations (|14[) are not manifestly hermi- 
tian; this can be remedied by introducing the variable 



VV = yjN-lifj a . 
Equations (| L4[) are then rewritten 



(15) 





% ~q^ ) ™ = -7V N - li^n - (1 + j)tp w + g\ip w \ 2 ip w ; (16a) 



d 



i-q^u, = - 1)% - jVN - lyj w + -JL-tytf^ 

(16b) 

with the initial state 



Mo))~Vn\ 1 )' 1 n 

For a large search space (JV 3> 1), Eqs. (fTB) reduce to 
d 

l Qi t l ) w ~ -jVNtpp - (1 + j)tp w + g\tp w \ 2 4> w ; (18a) 



*^Vfc « "7^m - 7^JV^ + jj\^\ 2 ^, (18b) 



with the initial state 



Mo) 
<M0) 



(19) 



Since ip^ and ip w are complex variables, they can be 
represented as 



/ N w e^, 



(20) 



where JV M and N w are the populations of bosons in the 
states and \w), respectively. Eqs. (ITS)) then corre- 
spond to four coupled nonlinear differential equations. 
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To reduce these to two coupled equations, one can define 
new variables 



r, = |^| 2 - |Vv| 2 = N w - N^; 



Eqs. (fT8"|) then become 



f] = 2jVN yji - rf sin(»; 



(21a) 
(21b) 

(22a) 



^ = 8- 9 -r l -2 1 ^N 7 J= cos(^) > (22b) 



where 8 is 



5 = 1 - Nj - 



(23) 



Eqs. (f2"2"j) are almost identical to those describing the 
Josephson dynamics of two weakly coupled Bose-Einstein 
condensates; see for example Eq. (2.6) in Ref. [36| . 
Note that taking the large- N limit is not necessary; the 
general- iV case is recovered simply by replacing N —> 
N — 1. The initial conditions for these variables are 

Ho),m) = (-^r- ) = ■ (24) 
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FIG. 1: The second set of fixed points as a function of g for 
TV = 1024, the dashed line is g+ and solid line is rj- and the 
thick black line is 770. 



1 . Fixed points 

The fixed points are obtained by setting Eqs. (|2"2"|) to 
zero. These are 



ii = 2m7r, raeZ; 770 = 0; 



(28) 



B. Complete search: 5 = 



and 



The complete search corresponds to the evolu- 
tion of the state from the initial state (p~7j) where 
{|^(0)| 2 , |^(0)| 2 } = corresponding to the 

maximum probability on the superposition of all un- 
marked sites, to the state with maximum probability on 
the marked site, i.e. {|^ M | 2 , \ipw\ 2 } = {jf^-Jf}- In 
other words, the search is complete when |r/"(t g )| = [f?(0)|. 
To determine if this is possible, it is convenient to inter- 
pret Eqs. (|22[) as Hamilton's equations of motion 



V 



8Hr 



8H, 



c 



drj 



(25) 



for some classical Hamiltonian He- It is straightforward 
to verify that a Hamiltonian satisfying both Eqs. (|25|) 
and J22]| is 



He 



8rj - |?7 2 + 2jVn^1 - if cos(0). 



(26) 



This has the same form of classical Josephson Hamilto- 
nian [33] • Since ^M- = 0, this Hamiltonian is a constant 
of motion. The value of He for a trajectory starting at 
the initial point {17(0), (j>(0)} = {— 1 + j^, 0} must there- 
fore be the same for the desired output {r](t s ), <f>(t s )} — 
{1 - f,0}. This is only possible if 28 (l - jf) =0 or 
(5 = 0, which corresponds to the homogeneous case. Set- 
ting 8 — for a complete search specifies a critical value 
for 7: 



7 



2-g 
2N 



(27) 



Note that the hopping coefficient must be positive, which 
requires g < 2. 



r?+ = -+ 
b 2 = (2m+ 1)tt; { 770 = 0; 



, _ 4( 3 -2)^ . 

a' 2 N > 



(29) 



_ 4(g-2 



As shown in Fig. [TJ r/ + and r\- approach zero as g de- 
creases and they both vanish at g = g* , where 
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2 + VN 



(30) 



To make further progress, one must identify the nature 
of the fixed points; a brief review of these concepts is 
given in Appendix [3] From Eqs. (|2"2"j) with (5 = 0, the 
functions appearing in the Jacobian (|A5[) are 



a(rj, (jy) = 2j*VNy / l-?f sin(</>); (31a) 

&(r?,^)=-|r7-2 7 *ViV J==cos(^). (31b) 
^ \Jl--q 1 

The Jacobian is therefore 

— r\ sin(0) (1 — ?? 2 ) cos(</>) \ 



J = 



2 7 ViV 



For the first set of fixed points (ry, 
Jacobian matrix becomes 







(0,2m7r) 



£ - £-£ 
2 VF 



(32) 

(0,27rar), the 



(33) 



5 




FIG. 3: Trajectory in phase space with N = 1024, g = 2g* 
and 7 = 7* = (2 — g*)/2N (light blue line). The red dots are 
?7_, the green ones are 77+, and the black ones are 770. The 
dark blue vectors are the eigenvectors of the Jacobian matrix 
at (0, (2m + 1)tt). 
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FIG. 2: Trajectory in phase space for TV = 1024, g = g* , 
7 = 7* = 2 ~E . The black and dark blue dots are marginally 
stable fixed points (g,4>) = (0,0) and (0, ±7r), respectively. 
Black and dark blue dashed circles are the orbits around these 
fixed points. The light blue line corresponds to the trajectory 
taken for a complete search. 



where the condition 7* 
The eigenvalues are 



(2 — g)/2N has been applied. 



A, = ± . f - g " 4 2)1 (34, 

Because TV » 1 and 7* > 0, the eigenvalues are strictly 
imaginary. The fixed points are therefore marginally sta- 
ble, or centers. The orbit in the vicinity of the fixed point 
at (77, (p) = (0, 0) is counterclockwise, as shown in black 
in Fig. dJ 

The second set of fixed points correspond to (77, (j>) — 
({?7±, 0}, (2m + 1)tt). Consider first the case (77, </5) = 
(0, (2m + 1)71"). The eigenvalues of the Jacobian matrix 
are then 



\±=±\ 



' (2-g)[-A + g(VN + 2)] 
2N 



(35) 



If S> 



these eigenvalues are both real. One is neg- 



2+vTV 

ative and the other positive, yielding an unstable saddle 
fixed point. The trajectories in the vicinity of these fixed 
points are depicted as dark blue vectors in Fig. [3] If 
g < 2+ 4 ^jy then both eigenvalues are imaginary, yield- 
ing a marginally stable fixed point or center. The orbits 
in the vicinity of the fixed points (77, </)) = (0, ±7r) are 
clockwise, as shown in black in Fig. [2] 

Consider next the fixed points (77, <j>) = (r]±, (2m+l)7r). 
For both cases, the eigenvalues of the Jacobian matrix are 



±- 



1 4(2- 9 )2 



TV 



(36) 



For g < 2J \j^ the solutions -q± do not exist. For 
g > — both eigenvalues become imaginary, yield- 
ing a marginally stable fixed point or center. The only 
difference between the two cases is that trajectories near 
the 77_|_ flow counterclockwise, opposite to the clockwise 
direction for those near r)- ; these are shown as green and 
red orbits in Fig. [3] 

Fig. [2] clearly shows that all of the fixed points in the 
g < g* regime are marginally stable or centers. Hence 
in this regime the trajectory starts from the initial point 
(77(0), 0(0)) = (— I + jt,0), rotates around the origin, 
and reaches the final point of the search (77 (t s ) , 4>(t a )) = 
(l — This behavior is depicted as the light blue 

line in Fig. [2] Thus, a complete search is attainable in 
this regime. 

In the other regime g > g* there is a saddle fixed point 
(77,0) = (?7_,7r) near the initial point (—1,0). The tra- 
jectory will continue along the positive </> and 77 will re- 
main close to 77-, so that it would never reach 7/ = +1. 
As shown in Fig. [31 it is not likely to have a complete 
search in this regime. In principle, one might still have 
a complete search for g > g* because the linearization 
procedure is strictly valid only right at the fixed points, 
but the range is not likely to be extensive. In any case, 
for any g £ [0, g*] a complete search is attainable. The 
value g — g* will be chosen for the remainder of the cal- 
culations. 



2. Search time 

Now that the nonlinear quantum search on the com- 
plete graph has been shown to be successful for a range of 
interaction strengths g, it is important to determine the 
scaling of the search time t s with the number of sites for 
large N. As shown in Fig.[5J the trajectory closely resem- 
bles a rectangle. Consider the dynamical equations (|2"2")l 
with 6 = and the large- N values of g = g* w 



6 



and 7 = 7* = (2 - g)/2N « l/N: 



V = — ^V 1 -V 2 sin(»; 



: COs(0) 



I]- 



(37a) 



(37b) 



For AT 3> 1, the right hand sides of Eq. (|3"Tf approach 
zero, and the trajectories are approximated by 77(f) ~ fci 
and ~ fc 2 , with /ci and /c 2 constants. These can be 
approximately decomposed into the following steps: 

(I) Constant 77: the initial point, (rj,(p) « (—1,0) — > 
(— 1, 0c), where C is the intersection of the trajec- 
tory with the axis; 

(II) Constant 0: (77, 0) w (-1,0 C ) (1,0c); 

(III) Constant 77: (77, 0) « (1,0 C ) ->■ (1,0). 

The search time f s can be written as 



t, 



dt+ dt 

(i) 7(ii) 



dt. 



(38) 



(in) 



For the first and third steps 77 is approximately constant, 
just as is approximately constant for the second step, 
so 



d4> 



1 dr] 
-1 V 



(39) 



The initial condition (pM|) gives a trajectory 77(f) s» — 1 
and — > 00 so that the first and third integrals make an 
insignificant contribution to t s . As shown in Fig. HI the 



77 m ki transition = 0—; 
w fc 2 transition 77 = — 1 - 
hovers in the vicinity of 0c] 
of a relaxation oscillator [36 

f s as 

_l 77 2 sin 
7T / 1 
2 \ sin(0, 



C is much faster than the 
1 (during which the phase 
This system is an example 
. Therefore one can express 



dry 



N 



)7_i 



A' 



(40) 



for large N. Besides the factor of l/sin(0 c ), this is the 
usual expression for the spatial search time. 

Setting 8 = 0, g = g*, and 7 = 7*, the classical Hamil- 
tonian He in Eq. (f2"o]l takes the form: 



2^/1 — rf cos(0) — ?7 2 



(41) 



For the initial condition rj(0) = — 1 + w — 1 for N 3> 
1, the classical Hamiltonian becomes -f/c ~ — anc ^ 
this value is preserved during the evolution. When the 



0.6 
0.4 
0.2 

-0.2 
-0.4 
-0.6 



20 



40 



100 , 120 




140 



FIG. 4: The phase as a function of time, when 5 = 0, 
S = S* = 7= and TV = 1024. 



trajectory depicted in Fig [5] crosses the axis at the 
point (?7,0) = (0,0c), the Hamiltonian is approximately 
He ~ 2cos(0 c )/v / iV. The value of the phase at this 
point is therefore C = cos _1 (— 1/2) = 27r/3. This is 
consistent with the time-evolution of the phase shown in 
Fig. |U The time for the nonlinear search in the large- N 
limit, Eq. (|40|) . is therefore 



ts = 77 



7T 2 
2 71 



7T 

N = ~F 



N, 



(42) 



slower than the linear search by a constant factor of 
2/V3 w 1.155. 

It is important to check that the linear search time 
is recovered in the case g — > 0. In this case one has 
7* = 1/AT (valid for all N). The trajectories will then 
cross the 77 axis at the initial condition (77(0), 0(0)) = 
(—1 + -h.O) when the classical Hamiltonian (l26l) takes 



2^1, valid for A^ > 1 



. When the tra- 
(0, C ), the clas- 



the value He 

jectory crosses the axis at (77, 0) 
sical Hamiltonian becomes He — —= cos(<fi c ). Because 
the Hamiltonian is a constant of the motion, one obtains 
0c ~ cos -1 (y^j ~ tt/2. The search time for the linear 

problem, Eq. (|4T)|) . is then t s = (n/2)\^N, consistent with 
expectations. 

C. Incomplete search: 5^0 



As discussed at the beginning of Sec. IIIIB1 a complete 
search is equivalent to finding a trajectory that makes the 
transition 77 = — 1 + — >• 1 — possible. Because He 
is a constant of the motion, however, a complete search 
requires 5 — 0. That said, it is conceivable that setting 
S =/= could yield a time t s where the relative occupation 
of the marked site would instead be 77 < 1 — . Consider 
again the fixed points of Eqs. (|22l) . Evidently = mix 
ensures that the right hand side of Eq. (|2^h 'l is zero, but 
77 = no longer accomplishes this for the right hand 
side of Eq. ([Hp) because 5 = 1 - N-/ - g/2 ^ by 
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assumption. Now that the search can not be complete, 
it is not necessary to keep the exact expressions for finite 
N, and one can work entirely in the limit N 3> 1. If 
one again sets 7 = 7* = 1/AT, then 6 = —g/2. Choosing 
g = a/N 13 , where a and (3 are both positive real numbers, 
gives S ~ which vanishes for large N. Eqs. (|22|) then 
become 



?7 = = — sin(0); 



AT 



aiX + rj) 
2N$ 



2 V 



: COs(</>). 



(43a) 



(43b) 



To determine the maximum value reached by 77, con- 
sider the classical Hamiltonian ([2"B|) which now has the 
form 



Mo.-SS&p + h^?^). (44, 



N 

At the initial conditon (rj,(p) ~ (—1,0), the classical 
Hamiltonian is Hc(0) = a/AN^ which is a constant of 
the motion. Note that for the incomplete search, the 
classical Hamiltonian is now positive. The accessible val- 
ues of (rj, (j>) are found by setting He = Hc(0). The 
A^-dependence disappears if /3 = 1/2 (i.e. g = a/y/W), 
and the only two real solutions correspond to 771 = —1 
and 



-1 



4(3x) 1 / 3 16 cos 2 (0) 
3a 2 {3X) 1 / 3 ' 



(45) 



where 



a 3 cos 2 (4>)V3y/27cP + 64 cos 2 (0) 



-I- 9a 4 cos 2 (</>). 



(46) 



Though the solution (|45| is a bit unwieldy, a few obser- 
vations can be immediately made. For small and large 
a, one obtains respectively 



mia < 1) 
m{u > 1) 



1 



■sec 2 ((/)); 



2cos 2 (0) 



1/3 



(47) 
(48) 



Note that only even powers of a appear in these expan- 
sions, indicating that the dynamics is unaffected by the 
sign of the nonlinearity. Eq. (|4"?| shows that a complete 
search is only possible for the non-interacting case a = 0, 
which is consistent with S — 0. For any finite interaction 
strength in this 5^0 regime, a complete search is not 
possible. The relative fraction on the marked site de- 
creases with a, and Eq. (|48|) reveals that it is asymptot- 
ically zero for very large nonlinearities (though one still 
requires g < 2). Physically, the large interaction between 
atoms favors the initial state which is the superposition 



of occupying all sites of the graph. This is the dynamical 
self-trapping which has been noted previously for inter- 
acting Bose- Einstein condensates [lH, H(| ■ 

Interestingly, rj^" 1 ^ is independent of the size of the 
search problem (keep in mind however that g — ot/\N 
so that the strength of the nonlinearity decreases with 
N). Note also that the maximum probability is reached 

for <p — 0. For example, ^j™ 3 *'' = 0.99 requires a s» 0.28. 
The case r?2 = can be obtained directly from Eq. (|45l) . 
and one obtains a = ±8cos(0). At this value of a, the 
probability of occupying the marked site exactly 1/2; for 
any smaller a it is larger. 

It remains to calculate the time for the incomplete 
search t s . As in the S — case, the right hand sides 
of Eqs. (|4"3")) approach zero, which means that the tra- 
jectories are approximated by constant lines. Again, the 
<f> evolution is much faster than the rj evolution for the 
r\ « rji — — 1 trajectory. The rj ~ r\i trajectory is not go- 
ing to be as fast because the \J\ — rfe term in Eq. (|43b ) 
is no longer almost zero. The search time is then approx- 
imately 



t. 



1 dr\ 
1 V 



d<j)VN 



2sin(^ c ) ' 7 0c _2M_ 2l 

2 



cos(0) 



(49) 



The critical angle 4> c is obtained by equating the classical 
Hamiltonian (gU) for (77, 0) = (O,0 C ) with H c (0); this 
gives a/4\/N = 2 cos((f> c ) / \/N or 



cos — 



7T 

2 



(50) 



for small a. For the linear case (a = 0), the critical angle 
coincides with that found in the previous section for the 
complete search. The first term in Eq. (|49p is therefore 
TrV^/2sm(4> c ) w (tt\/]V/2)/ v /1 - a 2 /64 which restricts 
the strength of the nonlinearity to |a| < 8 (recall that for 
I a I > 8 the probability on the marked state is less than 
one half). 

The second term can be simplified by assuming that 
the value of 77, given in Eq. (|45j) remains approximately 
constant over the range of integration {— (/> c ,0}. In fact, 
plotting 77 for a range of a shows that it varies from 
to its maximum value, Eq. (|47l) . but this increase occurs 
only over a very small region in the vicinity of (f> = (j) c . 
For small a, the time for the incomplete search, given by 
Eq. p9|) . can be approximated as 



N 1 + 



a 

128 



- VN 



d(i> 



1 



2tt Vl6 



a + 4cos($)/a 

(51) 



a 

128 



Because the a-dependent correction terms are negative 
in the regime < a < 8, the time for the incomplete 



8 



search is generally shorter than that for the complete 
search by a small factor dependent on the strength of 
the nonlinearity. 



IV. CONCLUSIONS 

In this work we considered the spatial search algorithm 
on lattice with the topology of the complete graph, under 
the assumption that the continuous-time quantum walk 
is effected by a zero-temperature Bose-Einstein conden- 
sate. In the mean-field approximation, the equations of 
motion become nonlinear and correspond to the discrete 
Gross-Pitaevskii equation. The analytical results, using 
methods in nonlinear dynamics, indicate that a complete 
spatial search remains possible even in the presence of 
nonlinearity. 

For a successful search, the nonlinear coupling con- 
stant must decrease with the system size as TV -1 / 2 and 
the inter-site hopping amplitude decreases as ./V -1 , where 
N is the number of sites. The latter condition coincides 
with the criterion for a complete search found for the lin- 
ear search problem on the complete graph [2(|. Under 
these circumstances, the search time is found to scale as 
t s oc y/~N, with an overall constant factor that depends 
weakly on the strength of the nonlinearity. The probabil- 
ity of success generically decreases with the strength of 
nonlinearity, but there are particular choices of param- 
eters where the success probability can be made unity. 
In the limit of zero nonlinearity, the present results com- 
pletely recover those of Ref . [20| . 

The dynamics of the nonlinear system agree closely 
with those of the linear quantum walk. This indicates 
that nonlinearity, as long as its strength is kept bounded 
for a given system size, is no impediment to the imple- 
mentation of a quantum spatical search. The results sug- 
gest that Bose-Einstein condensates consisting of huge 
numbers of particles can be candidates for the imple- 
mentation of useful quantum algorithms. The hopping 
amplitude and the strength of the nonlinearity need to 
be adjusted for a given size of the search space. In prac- 
tice, the former could be accomplished by adjusting the 
depth or spacing of the lattice [3jf . and the latter through 
the use of Feshbach resonances [30| • 

While the case of the complete graph addressed in 
this work is relatively straightforward to analyze the- 
oretically, it is not simple to produce experimentally. 
It would be preferable in practice to conduct the spa- 
tial search on a regular lattice, for example a square 
lattice in three or lower dimensions. Unfortunately, 
continuous-time quantum walks based on the ordinary 
discrete Schrodinger equation do not provide useful quan- 
tum speed-ups on these lattices over the classical search 
time, though discrete-time quantum walks can be con- 
structed that do [HI]. That said, if the particle dispersion 
relation is linear rather than quadratic, the full quantum 
spee d-up is achievable on a three-dimensional square lat- 
tice [21|. One practical strategy to achieve this is to 



artifically induce Dirac fermions by suitably preparing 
an optical lattice [4l|. More directly, the excitations of 
a zero-temperature weakly interacting BEC on a lattice 
are characterized by a linear dispersion relation, which is 
a hallmark of the underlying superfluidity in these sys- 
tems [33|. It is therefore conceivable that the nonlinear 
mean-field dynamics of the BEC on a cubic lattice would 
yield the full quantum speed-up for the continuous-time 
spatial search problem. This possibility will be explored 
in future work. 

During the preparation of this manuscript, we became 
aware of other work that investigated the same system 
addressed in the present study (42|. While their results 
are consistent with ours when there is overlap, in their 
work the strength of the nonlinear coupling is assumed 
to vary with time. Their methods and conclusions are 
therefore complementary to ours. 
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Appendix A: Characterizing fixed points 

This discussion follows Ref. [38|. Consider a general 
two-dimensional nonlinear system of equations: 

u = a(u,v); v — b(u,v). 

A fixed point or equilibrium point is defined by 

u 



(Al) 



Suppose (u* , v*) is a fixed point for this system, satisfying 
a(u* , v*) = b(u* , v*) = 0. Let u' = u — u* and v' = v — v*. 
For small u 1 and v' 

u' « a(u*,v*) + —{u*,v)u' + — (u*,v*)v'; A2a) 
ou ov 



v> » b(u*, v*) + ^(u*,v*)u' + ^-(u*,v*)v'. (A2b) 
ou ov 

Since (u*,v*) is a fixed point, a(u*,v*) = b(u*,v*) = 0. 
Close to the fixed points (u' < 1, v' < 1), Eqs. (|A"2l can 
be written as 



w' « Jw' , 



where 



w = 1 v' " 



(A3) 



(A4) 



9 



and J is the Jacobian matrix 

/ da_ d(L \ 

J = t M ■ (A5) 

The general solution of Eq. (|A3| when J is non- 
degenerate and invertible is 

w(t) = c lZl e Xlt + c 2 z 2 e X2 \ (A6) 

where \ and zi are the eigenvalues and eigenvectors of 
the Jacobian matrix </, respectively; Ci, C2 are constants 
which are determined by the initial conditions. There are 
four possibilities for the stability of the fixed points: 

1. Ai and A2 are both real: 

(a) Stable fixed point: If Ai < and A2 < 0, 

then w' — > as t — > 00. 

(b) Unstable fixed point: If Ai > and A2 > 0, 

then w' — y 00 as t — > 00. 

(c) Unstable saddle point: If Ai < < A2, 

then if w'(0) is a multiple of z\ then w' — > 



as t — y 00 (stable along direction of zi); al- 
ternatively if w'(0) is a multiple of z 2 then 
w' — > 00 as t — > 00 (unstable along direction 
of z 2 ). 

2. Marginally stable fixed point / center: Ai 

and A2 are both complex. If 5i(A;) = then 
\w'\ — > const as t — » 00. Trajectories circulate 
around the fixed point and eventually return to the 
initial point; these are closed orbits. 

The Hartman-Grobman theorem states that the dy- 
namics of the linearized system in the vicinity of hyper- 
bolic fixed points, where 7^ 0, will be similar to 
that of the original nonlinear system. If one or both 
eigenvalues don't satisfy this condition, the fixed point 
is non-hyperbolic and therefore the dynamics are frag- 
ile to the inclusion of nonlinearity. That said, suppose 
that (u*,v*) is an isolated fixed point that is marginally 
stable, and there exists a conserved quantity Hc{u,v). 
If (u*,v*) is a local minimum of Hc(u,v), then all the 
trajectories sufficiently close to (u*,v*) are closed. 
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